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The quantum self-trapping phenomenon of a Bose-Einstein condensate (BEC) represents a re- 
markable nonlinear effect of wide interest. By considering a purely dipolar BEC in a double-well 
potential, we study how the dipole orientation affects the ground state structure and the transition 
between self-trapping and Josephson oscillation in dynamics. Three-dimensional numerical results 
and an effective two-mode model demonstrate that the onset of self-trapping of a dipolar BEC can 
be radically modified by the dipole orientation. We also analyze the failure of the two-mode model 
in predicting the rate of Josephson oscillations. We hope that our results can motivate experimental 
work as well as future studies of self-trapping of ultracold dipolar gases in optical lattices. 

PACS numbers: 03.75.Lm, 34.20.Cf, 32.10.Dk, 32.80.Qk 



C3 



O 

> 

o> 
m 

00 

o 

> 

X 



I. INTRODUCTION 

Ultracold dipolar gases display a variety of unique 
properties that are absent in those dominated by S-wave 
scattering [l| . The successful Bose-Einstein condensation 
of 52 Cr atoms with a large magnetic dipole moment [2| 
and the first experimental realization of ultracold KRb 
polar molecular gas [3( are attracting even more experi- 
mental and theoretical interests in ultracold dipolar Bose- 
Einstein condensates (BECs) . Other than possible appli- 
cations in fields such as quantum information 0, Q , one 
main motivation along studies of ultracold dipolar gases 
is to explore new physics afforded by the anisotropic and 
long-range dipole-dipole interaction. For a recent review, 
see Ref. @. Indeed, by varying the shape of a dipo- 
lar BEC, the dipole polarization axis, and the trapping 
geometry 0, [H, Q , the partially attractive and partially 
repulsive dipole-dipole interaction of a dipolar BEC can 
be easily manipulated to a great extent. Even more dra- 
matically, by using the well-established Feshbach reso- 
nance technique [n| El El El, 

it is now possible to 
significantly reduce or completely shut off the effects of 
short-range interactions and hence realize purely dipolar 
gases, where the physics is completely dominated by the 
dipolar interaction. For example, superfluid transitions 
in purely dipolar Fermi gases [14 , 1 1 51 ] , fractional quantum 
Hall states in purely dipolar gas trapped in a rotating op- 
tical lattice , and collapse features of a purely dipolar 
BEC EH have been studied. 

In this work, we investigate the properties of a dipolar 
BEC confined in a double-well potential. We focus on 
the structure of the ground state wave function and the 
dynamical quantum self-trapping (QST) phenomenon 
which is among the most dramatic and counter-intuitive 
effects induced by the self-interaction of a BEC. The QST 



phenomenon has been extensively studied in non-dipolar 
BECs d Ei El and is closely connected with other 
physical contexts including the Josephson effect in super- 
conductors [11 and the superfluidity of 4 He [2l| ■ In par- 
ticular, the QST of a non-dipolar BEC in a double-well 
potential, which is sometimes called a Bose- Josephson 
junction, is now well understood fl8l I22I [H [24 1 and 
has been observed experimentally [25|. From a time- 
independent point of view, on the mean-field level the 
self-interaction of a BEC can induce the emergence of 
stationary states with a large population imbalance be- 
tween two wells. From a time-dependent point of view, 
the Josephson oscillation between the two wells can be 
suppressed by the self-interaction of a BEC, which can 
be regarded as a quantum Zcno effect because each atom 
in the BEC is being "measured" by all the other atoms. 

Our interest here is how the anisotropic nature of the 
dipolar interaction can be exploited to manipulate the 
dipolar BEC in a double-well potential, as manifested 
in the static properties of ground state structure and 
in the dynamical evolution of the system prepared out 
of equilibrium. We approach this problem by three- 
dimensional numerical simulations. We find a wide re- 
gion where the dynamics of a dipolar BEC can show 
either Josephson tunneling or QST without displaying 
signs of collapse. More significantly, we shall demon- 
strate that the dipole orientation can radically affect the 
transition from QST to Josephson oscillation. Further- 
more, to gain more insights into the dynamical QST, we 
construct an effective two-mode model which has been 
widely applied in non-dipolar systems. We find that this 
simple two-mode model successfully explains the numer- 
ically observed transition threshold between Josephson 
oscillation and QST, but fails to account for the oscilla- 
tion frequency in the Josephson oscillation regime. We 
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provide a thorough analysis of the success and the failure 
of the two-mode model. 

This paper is organized as follows. In Sec. II, after in- 
troducing a three-dimensional realistic model of a purely 
dipolar BEC trapped in a double- well potential, we study 
the structure of the ground state wave function. In Sec. 
Ill, we present detailed simulation results regarding how 
the dipole orientation affects the self-trapping of the sys- 
tem. In Sec. IV, we attempt to use a simple effective 
two-mode model to explain the dependence of the on- 
set of self-trapping of a dipolar BEC upon the dipole 
orientation. In Sec. V we discuss why our effective two- 
mode model cannot be used to estimate the oscillation 
frequency in the Josephson oscillation regime. Section 
VI concludes and summarizes this work. Appendix A 
presents some details regarding our numerical calcula- 
tions of the dipole-dipole interaction potential. 




II. GROUND STATE STRUCTURE 

To be specific, let us consider a dipolar BEC of 52 Cr, 
which has large magnetic moment /i = 6/ib (i^b is the 
Bohr magneton). We assume below the S-wave scatter- 
ing length is tuned to zero via Feshbach resonance, thus 
resulting in a purely dipolar atomic BEC. We further 
assume that the system is confined in a trapping poten- 
tial with a cylindrical symmetry in the x-y plane and a 
double- well structure along the z direction. That is, we 
consider a confining potential 

V(x, y, z) = y(^ 2 + cj 2 p y 2 + cj 2 z z 2 ) + Acxp(-z 2 /2a 2 ), 

where m is the atomic mass, to x = cu y = lu p is the confin- 
ing harmonic frequency in the transverse direction, uj z is 
the confining harmonic frequency in the longitudinal di- 
rection, A and Co are the height and the width of a Gaus- 
sian profile modeling a barrier between the two potential 
wells along the z direction. Without loss of generality, 
we assume that the dipole moments are polarized by the 
external magnetic field and are confined in the x-z plane. 
The two-body dipolar interaction potential is then given 
by 

U dd {r) = d 2 [r 2 - 2>{zcoblp + xsinip) 2 ] /r 5 , (1) 

where tp is the angle between the polarized dipole orien- 
tation and the z axis, d 2 — [i^ix 1 j^K with /io being the 
magnetic permeability of the vacuum. For convenience 
we define a dimensionless dipolar interaction parameter, 
D = (N — l)md 2 /(h 2 aho), where a^o = yS/ (muj z ) is the 
axial harmonic oscillator length and N is the total num- 
ber of atoms. We can then adopt a unit system where the 
units for length, time, and energy are given by ciho, 
and huj z , respectively. Unless specified otherwise, for the 
numerical results presented below we set uj p /uj z = 10, 
A = 4ftw 2 , (T = 0.2a/ lo , and D = 0.6 (which corresponds 
to about 2000 Cr atoms for the axial trapping frequency 
oj z = 2ir x 330 Hz). In such a confining potential, the 



FIG. 1: (Color online) Column density of the ground state 
density profile J ^(r)! 2 dx for two different angles (p. 
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FIG. 2: (Color online) The ground state energy, chemical 
potential as well as the population difference between the two 
wells as functions of tp. 

mean-field dynamics of a purely dipolar BEC at T = is 
described by the following dimensionless time-dependent 
Gross-Pitaevskii equation: 

i 9 -^ 1 = 4v 2 vKr,t) + V>)*( M ) 

+ J C/ (id (r-r')|*(r',t)| 2 vI/(r,i)dV,(2) 

where the macroscopic wave function ^ is normalized to 
unity. 

In our numerical simulations we adopt the time- 
splitting Fourier pseudospectral method [26|, Dis- 
cretization used in our calculations is At = 0.005 for 
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time, and Ax = Ay = 0.09375, Az = 0.25 in space. 
The calculation is performed in a box of size (L x , L yi 
L z ). The values of L x , L y , and L z are chosen such that 
the box is much larger than the size of the trapped con- 
densate. We impose zero boundary conditions for the 
wavefunction amplitudes at x — ±L x /2, y = ±L y /2, and 
z = ±L z /2. Typically we set L x = L y = 6, L z = 16. To 
evaluate the dipolar interaction potential that involves 
a convolution integral, we use the following fast Fourier 
transform technique, 

J U dd (r-r')\^(r\t)\ 2 d 3 r' = F-^FiU^F^ir,^ 2 }}, 

(3) 

where F and F 1 stand for fast Fourier transform and 
fast inverse Fourier transform, respectively. The term 
F[Udd( r )] m Eq. ([3]) is calculated analytically in the mo- 
mentum space (see detail in Appendix): 

F[U dd (r)] = J U dd {v)^d 3 v 

= nd 2 [2 sin 2 tp sin 2 9 k cos (2<^ k ) 

+ 2 sin(2iy9) sin(20k) cos 0k 

+ (4/3- 2 sin 2 95) (3 cos 2 k - l)] , (4) 

where (9 k and </> k are the polar angle and the azimuthal 
angle in the spherical coordinate system for the momen- 
tum space. 

The ground state is obtained by evolving Eq. flSj) in 
imaginary time. Figure [T] represents the ground state 
wave functions for two different dipolar angles: tp = 7r/36 
and tp = 57t/ 12. One can clearly see that in the former 
case, the wave function is asymmetric and localized in 
one of the wells; while for the latter, the wave function 
is symmetric and the population is equally distributed in 
both wells. Given the geometry of the trapping potential, 
the effective dipolar interaction is predominantly attrac- 
tive for small values of tp and repulsive for large values 
of tp. This explains the different structures of the ground 
state wave function at different angles. Figure [2] shows 
the ground state energy, chemical potential as well as the 
population difference between the two wells as functions 
of tp. It indicates that the symmetric/ asymmetric transi- 
tion as induced by the variation of tp is a continuous one 
and the critical angle is about tp — 0.22-7T. 

III. SELF-TRAPPING AND JOSEPHSON 
OSCILLATIONS 

We now turn to the dynamical properties of the sys- 
tem. We simulate the following situation. Initially we 
prepare the system in the ground state with tp = 7r/36. 
As shown in Fig. [TJa), for this dipole angle, the wave 
function is localized in one of the wells (for the example 
given, the atoms are localized in the left well). At t = 0, 
we suddenly change tp to some other value and study the 
ensuing dynamics of the system. 
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FIG. 3: (Color online) Column density of the condensate, (a) 
For the dipole orientation angle tp — tt/4, the population is 
mainly trapped in the left well instead of tunneling between 
the two wells, (b) Changing only the dipole orientation pa- 
rameter to tp = 57r/12, the tunneling between the two wells, or 
the Josephson oscillation, is observed. See the text for other 
system parameters. 
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In Fig. [3l we show the time evolution of V&(r) for two 
different final values of the dipole orientation parame- 
ter tp. In Fig. [3][a) tp is chosen to be 7r/4, and the 
condensate is found to remain localized in the left well 
for all times even though the ground state for p = ir/4 
should have population equally distributed in both wells 
[see Fig. [2jb)] . The system is hence clearly in the QST 
regime. By contrast, in the case of Fig. [3£b), tp is changed 
to 57r/12 and the population exhibits regular oscillations. 
The system in this case is clearly in the Josephson os- 
cillation regime. Because the only parameter difference 
between Fig. G2Ja) and Fig. E)Jb) is the final value of tp, 
Fig. [3] vividly demonstrates the tunability of the system 
through the orientation of the dipoles, a feature obviously 
absent in non-dipolar BEC's. 



To examine in more detail how the dipole orientation 
impacts on the QST, we next scan the value of tp and then 
plot the corresponding time evolution of a normalized 
population imbalance S between the two wells, namely, 
the population difference divided by the total number of 
atoms. Results are shown in Fig. |4j It can be seen that 
there exists a critical angle (p c f=a 0.377T: for (p < (p c , the 
system is in the QST regime; for tp > tp c , it enters the 
Josephson oscillation regime. We have also checked the 
critical regime in more detail by scanning tp in smaller 
steps. In particular, we also show in Fig. 0Ja) the popu- 
lation dynamics for tp — 13537r/3600, a value just above 
tp c . Interestingly, the oscillation dynamics for that case 
shows noteworthy critical oscillation behavior. Evidently 
then, our detailed results here further confirm that the 
dipole orientation of a dipolar BEC offers a simple and 
powerful means to manipulate the transition from QST 
to the Josephson oscillation. 



We have also studied many cases with other trapping- 
frequency ratios. For example, we let uj p /uj z vary from 
10.0 to 5.0 via increasing u> z . It is found that the fre- 
quency of the Josephson oscillation can be very sensi- 
tive to the trapping frequency ratio. For example, for 
lu p /lu z — 5.882, tp — 5n/l2, and the width of the initial 
wavefunction in z narrowed by -y/10/5.882, the Joseph- 
son oscillation frequency is more than 1.7 times larger 
than that shown in Fig. 2] (a) for the same value of tp. 
If Ljp/uj z — 5.882 and if the width of the initial wave- 
function in z is not narrowed , then with all other pa- 
rameters fixed, the condensate can already show clear 
signs of collapse (i.e., developing very high densities in 
our mean-field calculations) as we vary the dipole orien- 
tation. Interestingly, the condensate may collapse before 
we reach the critical point tp c for the transition between 
QST and the Josephson oscillation. This being the case, 
a stable dipolar BEC in a double-well potential may be 
always in the QST regime. 
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FIG. 4: (a) Population imbalance S(t) versus time, for dif- 
ferent values of the dipole orientation parameter tp. As tp 
exceeds a critical value tp c w 0.37-7T, the population imbalance 
starts to oscillate around zero, (b) Time averaged S(t) as a 
function of the dipole orientation parameter tp. When tp ex- 
ceeds tp c , the time averaged population imbalance suddenly 
decreases to zero. The transition between Josephson oscilla- 
tion and quantum self-trapping can hence be manipulated by 
tuning the dipole orientation. Other system parameters are 
the same as those used in Fig. 1. 



IV. A SIMPLE EFFECTIVE TWO-MODE 
MODEL 

To understand and gain more insights into the dynam- 
ical behavior presented in the previous section, we now 
construct a two-mode model that has been widely used 
for the study of QST for non-dipolar BEC's. To that 
end we assume that the main feature of the time-evolving 
wave-function \&(r, i) is captured by two normalized real 
basis states, 



tf(r,i) = ^i(t)$i(r) + Va(t)$2(r), 



(5) 



where ipi,2(t) = y/N^fie* 81 '*® , and $1,2(7) are local- 
ized in each of the two wells. The total number of atoms 
is given by \ipi | 2 + \ip2 1 2 = N T . Substituting Eq. © into 
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Eq. ([2]), we obtain 



eft 

•i [?/>i(i)V 2 $i(r)+V> 2 (t)V 2 $ 2 (r) 



+ [^i(*)^(r)*i(r)+^(t)V(r)$ 2 (r)] 
+ ||Vi| 2 Vi<fi(r) 1 [/ drf (r-r')|$i(r')| 2 dV 

+ |Vi| a lM*a(r) / C/ d(i (r-r')|$i(r')| 2 d 3 r' 

+ $i(r) /" (7 drf (r-r')$i(r')$ 2 (r')rf 3 r'] 

+^V 2 $i(r) /" U dd (r - r')$i(r')$ 2 (r')d 3 r' 

+ V-i|V- 2 | 2 [$ 2 (r) y C/ dd (r-r')$i(r')$2(r')tfV 

+ $i(r) y (7 drf (r-r')|$ 2 (r')| 2 dV] 

+ $^$ a (r) y C^,(r-r')$i(r')$a(r')d 3 r / 

+ |^ 2 | 2 ^ 2 $ 2 (r)y(7 dd (r-r')|$ 2 (r')| 2 d 3 rj. (6) 

As in previous studies of QST of non-dipolar BEC's, if we 
assume the overlap between the two modes being zero, 
namely, 

y ^(^(r^rwtfij; i,j = 1,2, 

J f (^(^(r^r^O; i^j; (7) 

then most of the terms in Eq. ([S]) will vanish. Indeed, 
such a great simplification is the main advantage of a 
two-mode picture in the first place. Adopting this zero- 
overlap assumption, we obtain two simple coupled equa- 
tions for ?/>i j2 (t): 

i^- = [(S0 + B 12 ) + ( J 4 11 -B 12 )|^i| 2 ]^i-^2, 

(8a) 

= [{E° 2 +B 21 ) + {A 22 -B 21 )\i; 2 \ 2 )]i !2 -KiP 1 . 

(8b) 

The parameters in the above two-mode equations of mo- 
tion are given by 

= yi[(V$ 4 ) 2 + l^(r)|V(r)]d 3 r, (9) 
A u = fu dd (v-r')\^ i (v)\ 2 \^ i (T , )\ 2 d 3 Td :i T', (10) 

Bi2 = B 2i= j C/ dd (r-r')|*i(r)| 2 |$ 2 (r')| 2 d 3 rdV, 

(11) 
(12) 



-(V$iV$ 2 ) + $iV(r)$ 2 



Here, Ajj and £?i 2 are the "on-site" and "off-site" contri- 
bution of the dipolar interaction term, respectively. The 
presence of the off-site interaction term B\2 is a unique 
feature for long-range interactions. Ef + B\ 2 can be re- 
garded as an effective on-site energy for the two wells, and 
An — B\ 2 can be regarded as an effective self-interaction 
strength, and k is the coupling strength between the two 
modes and hence gives the characteristic rate of oscilla- 
tions between the two modes if the nonlinear terms are 
neglected. Note that the expressions in Eq. ([5]) are very 
similar to those for a non-dipolar BEC in a double-well 
potential [l8| except that the off-site interaction term B\ 2 
is absent in non-dipolar systems. 

In terms of the population imbalance 5(t) = [| tjj-y (t) | 2 — 
\ r 4 > 2{t)\ 2 }/NT and a relative phase parameter <f>{t") = 
9-2.it) — 0i(t), Eqs. (|8a|) and (|8b|) assume an enlighten- 
ing form, 



S(t) = -Vl-S 2 (*)sh#(t)], 
0(t) = AE + AS(t) + Sl ' : ' 



(13a) 
cos[0(i)],(13b) 



where the time variable is rescaled by a factor 2k, with 
parameters AE and A defined by 

AE = ^l + (Au-A 22 )N^ 



A 



2k 4k 
(An + A 22 - 2B 12 )N T 
Ak 



(14b) 



Because Eqs. (I13ap and (|13b|) are exactly the same 
as those for a two-mode QST model of a non-dipolar 
BEC, it is clear that within the above zero-overlap ap- 
proximation, the underlying physics of the QST dynam- 
ics of a dipolar BEC should be similar to that of a 
non-dipolar BEC. Indeed, the mechanical analog of Eqs. 
(|13ap and (|13b[) is a classical non-rigid pendulum of tilt 
angle <fi, angular momentum 5, and a length proportional 
to yjl - S 2 {t) [18]. For an initial condition [5(0), <£(0)], 
this pendulum will oscillate around 5 = if and only if 
A < A c , where 



A c = 



2[1 - Aff5(0) + y/1 - 5(0) 2 cos 0(0)] 



(15) 



d J r. 



Translating back to our BEC context, we have that if 
A > A c , then the system will be in the QST regime; 
and if A < A c , then the Josephson oscillations can be 
expected. 

The above two-mode analysis also indicates that the 
QST phenomenon of a dipolar BEC can be easily ma- 
nipulated, because the important parameter A depends 
strongly on the dipole orientation as well as the trap ge- 
ometry. Furthermore, the expression of A in Eq. (|14b|) 
implies that the competition between the on-site interac- 
tion terms An and the off-site long-range term £?i 2 might 
play a role in the QST physics. For example, if under 
some circumstances the on-site contribution cancels out 
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the off-site contribution, then A will be small and hence 
QST will not occur. 

Our discussions so far are based on the strong as- 
sumption that a two-mode model still applies well to a 
purely dipolar BEC in a three-dimensional potential with 
a double- well structure. However, because our simula- 
tions show that the density profile in each well depends 
strongly on the population imbalance, a simple two-mode 
model is not expected to work satisfactorily. As such, we 
propose an effective two-mode mode constructed in a self- 
consistent manner, by first extracting useful information 
from our numerical simulations. Specifically, we use the 
following normalized Gaussian ansatz 

= r^— 3 , (16) 

to model the basis states 4>i j 2(r) in the two-mode model, 
with its parameters ai,bi,Ci,Xi,yi, and z% to be fitted by 
long-time average properties of the density profile in our 
full three-dimensional simulations of Eq. |2]). Once the 
parameters for the two basis states are obtained numer- 
ically, then the value of K can be obtained analytically, 
and the values of An and B\2, and hence the value of A, 
can all be obtained. In so doing, the location and width 
of the Gaussian ansatz for $j(r) may change with the 
initial condition, e.g., the initial population imbalance. 
This feature is hence outside the conventional two-mode 
model of QST. For that reason our self-consistent effec- 
tive two-mode model is expected to catch some features 
not available in a conventional two-mode model. 

The results from our effective two-mode model are pre- 
sented in Fig. [5] Figure |5£a) depicts how the values of 
the two-mode model parameters An, A 2 2, Bj.2 and k 
change with the dipole orientation parameter ip. Fig- 
ure [Sfb) shows the (^-dependence of A as well as A c . 
From Fig. GJb) it is seen that the value of A initially 
increases, and then decreases to almost zero at about 
ip f=a 0.377T. On the other hand, the value of A c , calcu- 
lated from Eq. (fT5|) with (f)(0) = (as an example), is 
seen to be less than A initially, and then exceeds A at 
about if w 0.397T. Therefore, as we tune the dipole ori- 
entation, our effective two-mode model displays a switch 
from A > A c to A < A c , thereby predicting the transition 
from the self-trapping regime to the Josephson oscilla- 
tion regime at a critical value ip ss 0.397T. The critical 
value obtained from the two-mode model matches with 
our previous numerical simulations of Eq. ([2]) , where we 
observed p c w 0.37-7T. We conclude that at least in our 
effective two-mode model, the anisotropic nature of QST 
for a dipolar gas can be understood in terms of the in- 
dependence of the two-mode model parameters. 

Results here indicate that when predicting the transi- 
tion point between the QST regime and the Josephson 
oscillation regime, the above-mentioned zero-overlap ap- 
proximation is acceptable. But it should be noted that 
this zero-overlap approximation is a very rough one, and 
there should be a high price for that. Indeed, as re- 
cently shown via a careful analysis of the validity of two- 
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FIG. 5: (Color online) (a) Values for the two-mode model 
parameters An, A22, S12, 

Ki , ELS functions of the dipolar ori- 
entation parameter ip, calculated by first fitting the Gaussian 
ansatz in (|16|) with the time-averaged properties of the den- 
sity profile in our three-dimensional simulations, (b) Value of 
A calculated from Eq. (|14b[l and value of A c calculated from 
Eq. (|15[) with (f>(0) — 0, as a function of the dipolar orienta- 
tion parameter ip. Note that for (p < <p c « 0.397T, A > A c ; 
and for ip > <p c , A < A c . 



mode approximations for one-dimensional systems [28|, 
the Josephson oscillation rate based on a simple two- 
mode theory may differ greatly from the actual rate. 
With that in mind we ask if our effective two-mode pic- 
ture can correctly describe other dynamical aspects. In 
particular, can our effective two-mode picture be em- 
ployed to predict a right order of magnitude of the fre- 
quency of the regular Josephson oscillations seen in the 
four bottom panels in Fig. HJa)? As it turns out, the 
answer is negative. For example, for <p = ir/2, the oscil- 
lation period observed from Fig. Ufa) is about t = 50; 
whereas the oscillation period from our effective two- 
mode model is about n/n ~ 0.157 [see Eqs. (|8a|) and 
(|8b[) ]. The difference is hence more than two orders of 
magnitude! Our next section is devoted to this interest- 
ing observation. 



V. DISCUSSIONS 

Here we provide a detailed analysis of the huge dis- 
crepancy between the two-mode parameter k and the ac- 
tual rate of the observed Josephson oscillations. Let us 
first re-examine the meaning of n defined in Eq. (|12[) . 
Clearly, for a three-dimensional problem, if the trans- 
verse direction is highly confined, the transverse motion 
will have a high kinetic energy and hence a large mag- 
nitude of wavefunction gradients with respect to x and 
y. According to Eq. (TT2"]) , this large transverse gradient 
will directly yield a large magnitude of k. Hence, for 
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tightly confined systems the magnitude of k is largely 
contributed by the kinetic energy in the transverse di- 
rection. In our two-mode treatment, k is identified as 
the parameter that determines the oscillation rate of the 
population between the two wells. As we show below, 
that k differs significantly from the actual oscillation rate 
obtained numerically originates from the zero-overlap ap- 
proximation made in the two-mode model. 

Specifically, if we now drop the previous zero-overlap 
assumption, Eq. (J6j> leads to 

«C 12 + £ff + B 12 (E° + B 21 )C 12 + k , 

l ~at = c£=i ^ + cl~i ^ 



y 12 



.dip 2 



kC 12 + E° + ggi 
C\ 2 - 1 



ty2 



(17a) 

{El + B 12 )C\ 2 + k . 

-Wi 



c\ 2 - 1 



(17b) 



where Q" and Q 2 are the nonlinear terms given by 



Qi 



f ' 



12 



1 



C?o - 1 



If - 



'12 " 

Ci 



^2 " 1 



o 2 



C 2 12 1 



o-i, 



°12 



-cr 2 



(18) 



with 



a! = (A n -S 12 )|Vi| 2 Vi+ 2£>i2^2 

+ I'ialV'ilVa + Gul^lVi 

<72 = {A 22 - B 12 )\lp2\ 2 1p2 + ID-iX^x 

+ D 2X ify\ + G 2 i\ipi\ 2 ip2 

+ GaiVa^ + ^ilV'ilVi, 



(19a) 



(19b) 



and 



C12 — 



$!(r)$ 2 (r)d 3 r, 



Aj - /"cT cW (r-r / )4 i (r)* i (r)|$ i (r / )| 2 d 3 rd 3 r', 
Gy = Ju M (T-i')$i(r)$ J (T)$ i (T f )$ J (i')d i TiPr', 
My = y"c/ dd (r-r')^(r)$,(r)|$ 3 (r')| 2 d 3 rdV. 

Evidently then, once we drop the zero-overlap assump- 
tion, the two-mode treatment can no longer give a simple 
picture of the dynamics. Nevertheless, Eqs. (|17ap and 
(|17b[) indicate that if we take into account the nonzero 
overlap C\ 2 , the effective coupling between the two modes 
will be quite different. Indeed, if we assume that in 
the Josephson oscillation regime the nonlinear terms Q" 1 
and Q 2 X have a rather small impact on the oscillation 
rate, we have the following corrected inter-mode coupling 
strength 



(£? + B 12 )C 12 + k 
G 2 , - 1 



(20) 



So long as CV,- is nonzero, this corrected coupling strength 
Kj can differ significantly from n. In particular, the ki- 
netic energy in the transverse direction will make a pos- 
itive contribution to E®Ci 2 but a negative contribution 
to K. As a result, these two terms will largely cancel each 
other and hence the kinetic energy in the transverse direc- 
tion will not directly enter into hi. For example, for the 
dipole orientation <p = 7r/2, using the same parameters 
as in our previous effective two-mode model, we obtain 
K\ ~ k 2 ps —2.5. One can then expect that the pre- 
dicted Josephson oscillation rate will be about one order 
of magnitude closer to our numerical value. Further tak- 
ing into account the nonlinear terms Qf l , one can expect 
the magnitude of the inter-mode coupling to be further 
reduced and hence better agreement with simulation re- 
sults can be obtained. For example, Q" 1 contributes a 
term — -Bi 2 |^2| 2 Gi2/(G 2 2 — 1) to the inter-mode coupling 
strength. This contribution will partially cancel the B\ 2 
term in K\ defined above. All these observations make 
it clear that the zero-overlap approximation is the main 
reason why our effective two- mode picture cannot be used 
to predict the rate of Josephson oscillations. 

Consistent with our simulation result that the Joseph- 
son oscillation rate depends strongly on the trapping ge- 
ometry, we also find that if we slightly change the width 
of the fitting mode wavefunctions $i j2 (r) and hence the 
value of C12, then the value of k may also change sig- 
nificantly. This implies that if we take Eqs. (|17ap and 
(|17b[) as an improved two-mode theory, then it is possi- 
ble to refit the mode wavefunctions and obtain the right 
Josephson oscillation rate. This approach is however not 
appealing to us, because by working with Eqs. (|17ap and 
(I17bp we lose the simplicity of a two-mode picture. 

So why our two-mode picture with the zero-overlap 
approximation can still give the correct transition point 
between the QST regime and the Josephson oscillation 
regime? The main reason lies in that our two-mode 
treatment already self-consistently used much informa- 
tion from the simulation results. In particular, in the 
QST regime, A c defined in Eq. (15) is dominated by AE 
and hence scales with k -1 , the same scaling behavior as 
A. As such, the actual magnitude of k will not affect 
the ratio of A/A c , the key index for a two-mode the- 
ory. In the Josephson oscillation regime, A c ~ 1 because 
AE ~ 0, and A << A c because of the large magnitude 
of k. This is also consistent with the condition for the 
Josephson oscillation regime in the two-mode theory. 



VI. CONCLUDING REMARKS 

To conclude, using three-dimensional numerical sim- 
ulations, we have investigated the structures of the 
ground state and the quantum self-trapping phenomenon 
of a purely dipolar Bose-Einstein condensate trapped 
in a double-well potential. The anisotropic nature of 
the dipole-dipole interaction is seen to impact dramat- 
ically on the transition between two dynamical regimes, 
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namely, the Josephson oscillation regime and the quan- 
tum self-trapping regime. This finding is the key result 
of this work. 

To gain useful insights we constructed a simple ef- 
fective two-mode model to understand the transition 
from the self-trapping regime to the Josephson oscilla- 
tion regime. Interestingly, though our two-mode picture 
is based on a rough zero-overlap approximation, the tran- 
sition point obtained from our two-mode picture is in 
excellent agreement with our numerical simulations. 

Somewhat expected, when it comes to the rate of the 
Josephson oscillations, our effective two-mode picture is 
no longer valid. We traced its failure to the zero-overlap 
approximation and discussed how we might be able to 
improve our effective two-mode treatment and find a 
better agreement with simulation results by lifting the 
zero-overlap approximation. However, by lifting the zero- 
overlap approximation, the simplicity of a two-mode pic- 
ture in understanding the quantum self-trapping is lost. 

There are now keen interests in the dynamics of BEC's 
in optical lattices. In particular, the quantum self- 
trapping effect is known to play an important role in 
determining whether or not a condensate can spread out 
in an optical lattice [29[ . The anisotropic nature of the 
quantum self-trapping of a dipolar BEC suggests that 
tuning the dipolc orientation can lead to the control of 
the self-trapping effect in optical lattices and hence the 
control of transport properties of dipolar BEC's in peri- 
odic potentials. We hope our work here may stimulate 
experimental efforts along this line. 
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(BX, WB and BL), National University of Singapore, 
by the start-up fund (WBS grant No. R-144-050-193- 
101/133) as well as the "YIA" fund (WBS grant No. 
R-144-000-195-123) (JG), National Univ. of Singapore, 
and by the NSF of US and the Welch Foundation (Grant 
No. C-1669) (HP). 



APPENDIX A: FOURIER TRANSFORM OF 
DIPOLE-DIPOLE INTERACTION POTENTIAL 

In this appendix we briefly outline how to use Fourier 
transformations to calculate the dipole-dipole interaction 
potential in our three-dimensional simulations. Assume 
that all the magnetic dipoles are aligned along an exter- 



nal magnetic field B(t), with 

B (t) = B [zcosip + sin ip (x cos a + y sin a)] , (Al) 

where x, y, z are units vectors in a cartesian coordinate 
system. Then the dipole-dipole interaction energy be- 
comes 



U dd (r) 



2 r 2 — 3 [z cos <p + sin tp (x cos a + y sin a) 



. ,—d 2 sin 2 (f cos 2 a 
-4v^F ^ (Yqo 



Y 20 ) 



/67r d 2 sin 2 <p sin(2a) 
U/ — ^ (Y 2 -2 - Y 2 2j 



6ir d 2 sin2</3COs(a) 




/6-7T d 2 sin 2ip sin(a) ,,, . 
m/^- ^—^(Y 2 - 1+ Y 21 ) 

Fk d 2 cos 2 ip d 2 sin 2 ip 

-4a / — ^ Y20 H ; > 

5 r 3 r a 



(A2) 



where Y\_ m is the standard spherical harmonics. Consider 
then the Fourier transform of U(r), 



U dd (k) = / U dd {v)e lkr d A v 



(A3) 



To evaluate this Fourier transform we first use the fol- 
lowing identity 



Z=0 m=-l 



where ji (x) is the spherical Bessel function, and 8^, 4>k 
are two spherical angles that define the direction of the 
k vector. Upon integrations over the radial coordinate r 
we finally have 

U dd (k) = ird 2 [2 sin 2 <p sin 2 6> k cos (20 k + 2a) + 
2 sin 2<p sin 28^ cos (0k + a) 
+ (4/3 - 2 sin 2 p) (3 cos 2 k - l)] . (A5) 

For our purpose here for which the external magnetic 
field is time-independent, we set a = and hence obtain 
Eq. ®. 
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